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ABSTRACT 

We investigate a numerical solution of the flavor-nonsinglet Altarelli-Parisi equation 
with next-to-leading-order a s corrections by using Laguerre polynomials. Expanding 
a structure function (or a quark distribution) and a splitting function by the Laguerre 
polynomials, we reduce an integrodifferential equation to a summation of finite number 
of Laguerre coefficients. We provide a FORTRAN program for Q 2 evolution of non- 
singlet structure functions (Fi, F 2 , and F 3 ) and nonsinglet quark distributions. This 
is a very effective program with typical running time of a few seconds on SUN-IPX or 
on VAX-4000/500. Accurate evolution results are obtained by taking approximately 
twenty Laguerre polynomials. 
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Program Summary 



Title of program: LAG2NS 

Computer: SUN-IPX (VAX-4000/500); Installation: Saga University Computer Center 
(The Research Center for Nuclear Physics in Osaka) 

Operating system: SUN-OS-4.1.2 (VAX/VMS V5.5-2) 

Programming language used: FORTRAN 77 

Peripherals used: laser printer 

No. of lines in distributed program, including test data, etc.: 1111 

Keywords: Altarelli-Parisi equation, numerical solution, Q 2 evolution, Laguerre poly- 
nomials. 

Nature of physical problem 

This program solves the Altarelli-Parisi equation for a spin-independent flavor- nonsinglet 
structure function or quark distribution. 

Method of solution 

We expand an initial quark distribution (or a structure function) and a splitting func- 
tion by Laguerre polynomials. Then, the solution of the Altarelli-Parisi equation is 
expressed in terms of the Laguerre expansion coefficients and the Laguerre polynomi- 
als. 

Restrictions of the program 

This program is used for calculating Q 2 evolution of a spin-independent flavor-nonsinglet 
structure function or quark distribution in the leading order or in the next-to-leading- 
order of a s . The double precision arithmetic is used. The renormalization scheme is 
the modified minimal subtraction scheme (MS). A user provides the initial structure 
function or the quark distribution as a subroutine. Examples are explained in sections 
3.2, 4.13, and 4.14. Then, the user inputs thirteen parameters explained in section 3.1. 

Typical running time 

Approximately five (three) seconds on SUN-IPX (VAX-4000/500) if the initial distri- 
bution is provided in the form of aix 6l (l — x) Cl + a2X b2 (l — x) c ' 2 + • • • in the subroutine 
GETFQN. If Laguerre coefficients of the initial distribution in the subroutine FQNS(x) 
are calculated by a GAUSS quadrature, the running time becomes longer depending 
on the function form. 



2 



LONG WRITE-UP 



1. Introduction 



Structure functions measured in deep-inelastic lepton-nucleon scattering depend in 
general on two kinematical variables, Q 2 = —q 2 and x = Q 2 /2P ■ q, where q is the 
four-momentum transfer and P is the nucleon momentum. If the nucleon is pointlike, 
the structure functions depend only on the variable x. This assumption is called 
Bjorken scaling hypothesis. However, observed experimental data indicate weak Q 2 
dependence, and this fact is referred to as scaling violation. Although the structure 
functions themselves could not be calculated exactly in QCD except for lattice QCD 
methods, the scaling violation can be evaluated in perturbative QCD. The scaling 
violations of the structure functions have been investigated extensively, and an intuitive 
way to describe the phenomena is to use the Altarelli-Parisi equation pj: 
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where q NS (x,Q 2 ) is a flavor-nonsinglet quark distribution, a s (Q 2 ) is the running cou- 
pling constant in Appendix A, and P^s(x) is the nonsinglet splitting function in Ap- 
pendix B. The above equation describes the process that a quark with the nucleon's 
momentum fraction y radiates a gluon and it becomes a quark with the momentum frac- 
tion x. The splitting function Pns(z) determines the probability for a quark radiating 
a gluon such that the quark momentum is reduced by the fraction z. Next-to- leading- 
order (NLO) corrections to the Altarelli-Parisi equation are included in the coupling 
constant a s and in the splitting function P NS 
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q NS (y,Q 2 ), (1.2) 



where the running coupling constant is the NLO expression af LO in Eq. (A. 2), and 
the leading order (LO) and NLO splitting functions are given in Eqs. (B.l), (B.3), 
and (B.5). The leading order and the next-to-leading order are abbreviated to LO and 
NLO respectively throughout this paper. 

Because the above Q 2 evolution equation is very important for testing the pertur- 
bative QCD and is often used theoretically and experimentally, it is worth while having 
a computer program of solving it accurately without consuming much computing time. 
There are several methods for solving the above integrodifferential equation. These 
include brute-force methods ||, Mellin-transformation methods ||, methods of using 
orthogonal polynomials 0|, and others. Among these, the Laguerre method in Refs. 
[[|, |5|, |6| is considered to be very effective. In Ref. 0, we provide the Laguerre method 
program LAG1, which deals with the Q 2 evolution in the leading order. The purpose 
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of this paper is to extend the nonsinglet program in Ref. || by including the NLO 
corrections. 

In section 2, we explain the Laguerre method for solving the evolution equation 
with the NLO corrections. In section 3, we supply necessary information for running 
our program LAG2NS. The details of each subroutine are discussed in section 4, and 
numerical results are explained in section 5. Explicit expressions for the running cou- 
pling constant, the splitting function, the Laguerre polynomials, and other necessary 
functions are given in appendices. 



2. Laguerre polynomial method 

The Laguerre polynomial method was first investigated by Furmanski and Petronzio 
[|J. The quark distribution and the splitting function are expanded in the Laguerre 
polynomials, then the solution is given by the Laguerre expansion coefficients and the 
Laguerre polynomials. Numerical analysis of the method was studied by Ramsey |J . 
Especially, the fast-computing Laguerre method was developed if the initial distribu- 
tion is expressed in a simple analytical form ax b (l — x) c . According to Ref. f5[, the 
Laguerre polynomial method is considered to be more effective than other methods 
in convergence and in computing time. This is because the first several polynomi- 
als resemble parton distributions. Based on these investigations, a useful FORTRAN 
program was published in Ref. ||. However, the program can be used only for the 
evolution in the leading order of a s . We discuss the solution for the evolution equation 
including the NLO corrections in the following. We refer the reader to the papers in 
Refs. ||, ||, H for more complete account of the Laguerre method. 

We outline the Laguerre method with the NLO corrections. In the following, we 
use the NLO running coupling constant unless we specify LO. We use MS as the 
renormalization scheme. Because the quark distribution multiplied by x satisfies the 
same evolution equation, we define q(x) = xq(x) and P(x) = xP(x) and rewrite Eq. 
(1.2) as 
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In order to remove the Q 2 dependence in front of the integral, we use the variable t 
defined by 
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instead of Q 2 , where the running coupling constant a s is given in Eq. (A. 2) and the 
constant (3q is in Eq. (A. 3). For example, Q 2 evolutions from Qq=4 GeV 2 to Q 2 =10, 
100, and 1000 GeV 2 with A = 0.2 GeV and four flavors correspond to t =0.00384, 
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0.112, and 0.166. Then, the integrodifferential equation becomes simpler 
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where R(x) is given by P(°)(x) and 
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The constant Pi is given in Eq. (A. 3). 

We define the evolution operator E(x,t) by the convolution form 



1 ^-E(*,t)q NS (y,t = 0) 
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(2.5) 



Namely, the function E(x,t) describes the Q 2 evolution from the initial distribution 
at Qq (t = 0) to the distribution at Q 2 . From Eqs. (2.3) and (2.5), it satisfies the 
integrodifferential equation 
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Because Eq. (2.6) and Eq. (2.3) are exactly the same form, it may seem unnecessary 
to define such operator. The introduction of the evolution operator simplifies the 
following formalism because it is the S function at t — 0, E(x,t = 0) = 6(1 — x). 
Laguerre coefficients of E(x,t = 0) are calculated by using Eq. (C.3), and they are 
E n {t = 0) = 1 for all n. We write the function E(x,t) by LO and NLO contributions 
as 

(2.7) 
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We substitute this equation into Eq. (2.6). Then, considering that the LO evolution 
operator E^(x,t) satisfies the LO Altarelli-Parisi equation 
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we find that the NLO operator satisfies the equation 
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£ (1) (x,t) 



Pa 



1 - 



a s (t) 
a s (0) 



f d -lR( X -) E^y,t) 
Jx y \yj 



(2.10) 



5 



So we do not have to solve equations for both E^°'(x,t) and E^'(x,t). We solve Eq. 
(2.8) for getting E^°'(x,t), which is then substituted into Eq. (2.10) for calculating 
EV>(x,t). 

We introduce variables x' and y' by x' = — In x and y' = — In where < x, y < 1, 
because the Laguerre polynomials are defined in the region < x' < oo. Using 
these variables, we expand the functions E(x = e~ x ) and P(x = e~ x ) in terms of 
the Laguerre polynomials, f(e~ x ) = f n L n (x'). Then our problem is to obtain the 

n 

expansion coefficients E n in terms of the coefficients of the splitting function P n . The 
details of solving Eq. (2.8) are discussed in Ref. || as well as in Refs. [|, and we 
summarize the solution in Appendix C. The NLO evolution operator is then given by 
using the calculated E^(t) 
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Ri is the Laguerre coefficients of R(x). R^i is defined as 
operators in Eqs. (C.5) and (2.11), we finally obtain 
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In this way, we reduced the original integrodif- 



ferential equation to a sum of finite number of Laguerre expansion coefficients and the 
Laguerre polynomials. If the Laguerre coefficients of the splitting functions and the ini- 
tial distribution are calculated fast, the Altarelli-Parisi equation is solved numerically 
without consuming much computing time. 

We have discussed a solution for the Q 2 evolution of a nonsinglet quark distribu- 
tion. Q 2 evolution of a nonsinglet structure function is calculated in a similar way. 
The only modification is to take into account NLO corrections from the coefficient 
function. A nonsinglet structure function Fns{x,Q 2 ) is expressed as a convolution of 
the corresponding nonsinglet quark distribution q NS (x, Q 2 ) and the coefficient function 
C NS (x, a s ) as §, U 



F NS (x,Q 2 ) = [ — C NS {-,a 
Jx y y 

The coefficient function is 

C NS (x,a s ) = 5(1 
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B s (x) is the NLO correction, and those for the structure functions (F±, F2/X, and F 3 ) 
are listed in Eq. (B.6). Combining the evolution equation for the quark distribution 
in Eq. (1.2) with the convolution in Eq. (2.14), we obtain the Q 2 evolution equation 
for the nonsinglet structure function as M 
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This is exactly the same equation with Eq. (1.2) if we replace Pj^l(x) by P- 
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20 o B NS (x). Therefore, we do not have to solve the above equation for F N s(x,Q 2 ) 
independently. The solution Eq. (2.13) can be used with the replacements 

q NS (x, Q 2 ) -> xFx(x, Q 2 ), F 2 (x,Q 2 ), or xF 3 (x,Q 2 ) , (2.17a) 
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In the following sections, we explain our program of calculating Eq. (2.13) or the 
corresponding equation for the structure function. 
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3 Description of input parameters 
and input distribution 



3. 1 Input parameters 



There are thirteen input parameters in the main program. 



(1) 
(2) 



(5) 
(6) 
(7) 
(8) 
(9) 



FO 

IORDER 



number of quark flavors (F0=3 or 4). 



(3) ITYPE 



(4) IMORP 



1, 
2, 
1, 
2, 
3, 
4, 

1, 
2. 



leading order in a s ; 
next-to-leading order, 
structure function xF 1 



NS , 



2 

jNS , 



Q02 
Q2 

DLAM 
NMAX 
IQN 



'(x,Q 2 ); 
xF-(x,Q 2 ); 
quark distribution xq NS (x, 
q — q type distribution; 
q + q type distribution, 
initial Q 2 (= Qq) at which an initial distribution is supplied. 
Q 2 (in GeV 2 ) to which the distribution is evolved. 
QCD scale parameter Aqcd in GeV. 

maximum order of Laguerre polynomials (0<NMAX<31). 



= 1, calculate the Laguerre coefficients of the initial distribution by 
Eq. (C.4); 

= 2, calculate them by a Gauss-Legendre quadrature. 
= 1, linear scale in x; 
= 2, logarithmic scale in x. 
XMIN = minimum of a; (x ^ 0). 
XMAX = maximum of x. 

NXSTEP = number of x at which distributions are calculated. 



(10) ILOG 

(11) 
(12) 
(13) 

The above ILOG, XMIN, XMAX, and NXSTEP are parameters for writing evolu- 
tion results as the function of x. For example, the evolved distributions are calculated at 
x=0.1, 0.2, 1.0 for ILOG=l, XMIN=0.1, XMAX=1.0, and NXSTEP=9. XMIN=0 
should be avoided due to a numerical problem. In the end of this paper, a test run result 
is given for the parameter set: F0=4.0, IORDER=2, ITYPE=4, IMORP=l, Q02=4.0, 
Q2=200.0, DLAM=0.19, NMAX=30, IQN=1, ILOG=2, XMIN=0.01, XMAX=1.0, 
and NXSTEP=30. 



3.2 Input distribution 

If the initial distribution is written as a\x 1 (l-x) ci +a 2 x b2 (l-x) C2 + ••■, the input 
IQN=1 should be chosen in order to save computing time. In this case, the subroutine 
GETFQN should be supplied to calculate Laguerre coeficients of the input distribution 
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at Qq, [xq NS (x)] n or [xFNs{x)] n . In this subroutine, the nth Laguerre coefficient for 
the initial distribution ax b (l — x) c is calculated by calling the function QN(n,a,b,c). 
An example is given in Sec. 4.14. 

If the initial distribution cannot be written in the simple form, the input IQN=2 
should be chosen. Then, the initial distribution at Ql should be supplied as the double 
precision function FQNS(:r). An example is given in Sec. 4.13. 

4. Description of the program LAG2NS 

4 . 1 Main program LA G2NS 

The main program reads thirteen input parameters in section 3.1 from the input 
file 5. t defined in Eq. (2.2) is calculated in LO and in NLO. Then, the subroutine 
GETFQNS is called for calculating the Q 2 evolution. 

4.2 Functions ALPHA0(Q2) and ALPHA 1(Q2) 

Running coupling constants in LO and NLO are given in ALPHAO and ALPHA1 
respectively. The NLO a s is given in the MS scheme. See Appendix A for the details. 

4.3 Subroutines GETFQNS, OBTAIN 1, OBTAIN2, and OBTAINF 

The subroutine GETFQNS calculates Q 2 evolution of a nonsinglet distribution. 
First, the subroutine OBTAIN1 or OBTAIN2 (OBTAINF in the structure function 
case) is called, depending on the input parameter IORDER. These subroutines call all 
special functions used in the program. If IORDER=l, the subroutine OBTAIN1 calls 
necessary functions for calculating the LO evolution. Then, the evolved distribution 
in the leading order is calculated by calling the function FQNSO(x) for Bjorken x 
from XMIN to XMAX. The evolved distribution is written on the output file 6. If 
IORDER=2, OBTAIN2 (or OBTAINF) is called and the evolved distribution including 
the NLO corrections is calculated in the same way by calling the function FQNSl(x, t). 
The subroutine OBTAIN2 (OBTAINF) calls necessary functions for calculating NLO 
evolution results. 

44 Subroutines GETFACT, GETZETA, and GETZFUN 

These are necessary functions for calculating the LO evolution. See Ref. || for the 
details. 

4.5 Subroutine GETPNO 

This subroutine calculates Laguerre coefficients of the splitting function in the lead- 
ing order; coefficients are stored in the array PNSNN(n + l)[n = — NMAX] = 
(xPjsis)n- Then P n —P n _i is calculated and results are stored in the array PTILNSN (n+ 
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1) [n = O-NMAX] . These are used for calculating B k n in the subroutine GETBKM1. 

4.6 Subroutines GETRN3M, GETRN4M, GETRN3P, and GETRN4P. 

These subroutines store Laguerre coefficients of the splitting function R(x) [Eq. 
(2.4)] in the array RN(n + 1) [n = — NMAX]. These coefficients are numerically 
calculated in a separate program. GETRN3M is for evolution of a u q — q" type distribu- 
tion in the three-flavor case, and GETRN4M is in the four-flavor case. GETRN3P and 
GETRN4P are for a u q + q" type distribution with three and four flavors respectively. 
These are used for calculating E$ [Eq. (2.12)] in the subroutine GETEMT1. See Ref. 
H for the explanation of the u q — q" and u q + q" type distributions. 

4.7 Subroutines GETRN and GETFN 

In these subroutines, R n — R n -i is calculated and results are stored in the array 
RTILN (n + 1) [n = - NMAX] . These are used for calculating Eg) [Eq. (2.12)] in 
the subroutine GETEMT1. In the case of structure-function evolution, R(x) includes 
the coefficient function B NS (x) as shown in Eq. (2.17b). R(x) defined in Eq. (2.4) is 
used in GETRN, and Eq. (2.17b) is used in GETFN for calculating R n — R n -\. 

4.8 Subroutines GETF1N, GETF2N, and GETF3N 

Laguerre coefficients of the function B NS (x) are supplied in these subroutines. 
These coefficients are numerically calculated in a separate program. The functions 
B NS (x) for the structure functions, Fx, F 2 /x, and F 3 , are given in Eq. (B.6) and 
their (actually B NS (x) multiplied by x) Laguerre coefficients are supplied in GETF1N, 
GETF2N, and GETF3N respectively. The Laguerre coefficients are stored in the array 
FNN (n + 1) [n = - NMAX] . 

4.9 Subroutine GETBKM1 

This subroutine calculates in Eqs. (C.6) and (C.7), then they are used for 
calculating the LO evolution operator in Eq. (C.5). 

4-10 Function QN(N,a,b,c) and subroutine GETQND 

QN(N,a,b,c) is used for calculating Laguerre coefficients of the initial distribution 
ax b {l — x) c in the subroutine GETFQN in the end of this program. 

GETQND calculates Laguerre coefficients by using the Gauss quadrature subrou- 
tine DGAUSS with the accuracy of 1CT 8 . 
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Ill Subroutines GETEMTO and GETEMT1 



GETEMTO calculates Laguerre coefficients of the evolution operator (t) in Eq. 
(C.5). In calculating this equation, we need [xPns]o calculated in the subroutine 
GETPNO and B k m in the subroutine GETBKM1. 

To calculate the evolution of a quark distribution, the differences E^\t) — E^\{t) 
are calculated; the resulting differences are stored in the array EPSNS(n + l)[n = 
0-NMAX). 

GETEMT1 calculates Laguerre coefficients of the evolution operator E$p(t) in Eq. 
(2.12). To calculate this, RTILN calculated in subroutine GETRN (GETFN) and 
EPSNS are required. Then, the differences E^\t) — E^l^t) are calculated and the 
results are stored in the array EPSNSl(n + l)[n = — NMAX}. They are used for 
calculating the evolution of a quark distribution in the next-to-leading order. 

4.12 Functions FQNSO(x) and FQNSl(x) 

The evolved distribution in the LO (in the NLO) is calculated by using Eq. (2.13) in 
the function FQNSO (FQNS1). It is the quark distribution FQNSO = xq NS (x,Q 2 ) for 
the input parameter ITYPE=4, and it is the structure function FQNSO = Fns(x, Q 2 ) 
[xFf s , F^ s , or xFf s ] for ITYPE=1, 2, or 3. 

4.13 Function FQNS(x) 

If the initial distribution cannot be expressed in the simple form ax b (l — x) c 
(IQN=2), the distribution should be supplied as the function FQNS(x). The HMRS-B 
valence-quark distribution xu v (x) + xd v (x) |TtJ and the CCFR F 3 structure function 
xF 3 (x) |Tl| are given as examples 

[xu v + xd v ] HMRS ~ B = 0.5469x a237 (l - x) 4 - 07 (l + 23.8x) , (4.1) 

[xF 3 } CCFR = 5.976x°' 766 (l - x) 3 - 101 . (4.2) 



4-14 Subroutine GETFQN: This subroutine or the function FQNS(x) should be sup- 
plied by a user. 

If the initial distribution can be expressed in the form aia; fel (l — x) Cl + a2X b2 (l — 
x) C2 + ■ ■ • (IQN=1), the Laguerre coefficients are calculated by using the function 
QN(N,A,B,C) as QN(N, a u b u a) + QN(N, a 2 , 6 2 , c 2 ) + ■ ■ ■. The HMRS-B and CCFR 
distributions are given as examples 

[xu v + xd v )* MRS - B = QN(n, 0.5469, 0.237, 4.07) + QN(n, 13.016, 1.237, 4.07) , (4.3) 

[xF 3 }^ CFR = QN(n, 5.976, 0.766, 3.101) . (4.4) 



11 



5. Numerical analysis 



Our FORTRAN-77 program LAG2NS can be run in double precision arithmetic. 
This is a program for the Q 2 evolution of a spin-independent flavor-nonsinglet structure 
function or quark distribution with or without the NLO corrections. 

The initial distribution or its Laguerre coefficients should be provided by the func- 
tion FQNS(x) or the subroutine GETFQN within the program. If the initial distribu- 
tion can be expressed in the analytical form a\X 1 a user 
selects the input IQN=1 and provides the subroutine GETFQN. In this case, comput- 
ing time is very short. However, if the distribution cannot be expressed in the simple 
analytical form, the function subroutine FQNS(x) should be supplied. Our program 
calculates Laguerre coefficients of this distribution by a GAUSS quadrature with the 
accuracy of 10~ 8 . So it could take a significant amount of computing time depending 
on complexity of the initial function. Then, thirteen input parameters should be pro- 
vided from the input file 5 for running the program. Evolution results are written in 
the output file 6. An example of the output is shown in the end of this paper as TEST 
RUN OUTPUT. 

We tested our program in comparison with other publications. In order to save 
computing time, the Laguerre coefficients of the NLO splitting functions and the co- 
efficient functions are calculated in separate programs. Obtained Laguerre coefficients 
are supplied in our program as discussed in sections 4.6 and 4.8. In the separate pro- 
grams, the splitting-function and coefficient-function subroutines are tested by calcu- 
lating their moments, namely anomalous dimensions. Calculated results are compared 
with those in Refs. [12|, Because the NLO expressions in the appendix B are 
rather complicated, the Laguerre coefficients are evaluated numerically by the Gauss 
quadrature with the accuracy of 10~ 12 . We compared our NLO evolution results of 
u v (x,Q 2 ) + d v (x,Q 2 ) with those of the HMRS-B program |nj. The HMRS-B distri- 



bution at Q 2 =4 GeV 2 is evolved to Q 2 =200 GeV 2 with four flavors and 
GeV. We find that differences between our evolution and the HMRS one are of the 
order of 0.5% in the x range 0.01 < x < 0.5. The difference becomes slightly larger in 
the large x region; however, it is still of the order of 2%. Hence, our program results 
essentially agree with the HMRS Q 2 evolution. We also checked our NLO corrections 
in comparison with Ref. |14| . Namely, we choose the starting nonsinglet F 2 structure 
function F 2 (x,Q 2 = 5GeV 2 ) = 3.8075x°- 56 (l - x) 2 - 71 . Then the structure function at 
Q 2 = 200 GeV 2 is calculated with A=0.3 GeV and four flavors. The results agree with 
those in Fig. lb of Ref. jnj; however, we find a slight deviation in the medium x 



region. 

Next, we show examples of our evolution results. We choose the HMRS-B dis- 
tribution at Q 2 =4 GeV 2 [u v + d v = 0.5469x a237 (l - x) 4 - 07 (l + 23.8a:)] as the initial 
distribution. This distribution is evolved to the one at Q 2 =200 GeV 2 by assigning the 
input parameters F0=4.0, IORDER=2, ITYPE=4, IMORP=l, Q02=4.0, Q2=200.0, 
DLAM=0.19, NMAX=10, 20, or 30, IQN=1, ILOG=2, XMIN=0.01, XMAX=1.0, and 
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NSTEP=50 and by providing the subroutine GETFQN. Evolved valence-quark distri- 
butions are shown in Fig. 1. We find that reasonably accurate results are obtained by 
taking about twenty Laguerre polynomials. The accuracy is within 3% in the x range 
0.001 < x < 0.8. However, there is a tendency to become worse in the very small x 
region (x < 0.001) and in the very large x region (x > 0.9). 

We compare our evolution results with the CDHSW neutrino data in Fig. 2. 
We choose the CCFR F 3 (x) distribution at Q 2 =3 GeV 2 as the initial distribution in 
both LO and NLO cases in order to see NLO contributions. It is given by xF 3 (x, Q 2 = 
3 GeV 2 ) = 5.976x°' 766 (l - x) 3 - 101 . We calculate Q 2 variations with A=0.21 GeV, four 
flavors, and the q — q splitting function by providing the subroutine GETFQN. The 
Q 2 variations are calculated at x=0.045, 0.225, and 0.55. The NLO contributions are 
important at small Q 2 (~ 1 GeV 2 ) and at small x as shown in Fig. 2. Our Q 2 variations 
are consistent with the CDHSW experimental data. 

Typical running time is a few seconds on SUN-IPX or on VAX-4000/500. Consid- 
ering the accuracy and the short running time, we conclude that the Laguerre method 
is considered to be very effective for numerical solution of the Altarelli-Parisi equation. 

We have investigated a numerical solution of spin-independent flavor-nonsinglet 
Altarelli-Parisi equation. We will work on flavor- singlet structure functions and also 
on nuclear structure functions. 
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Appendix A. Running coupling constant 



We give explicit expressions of the running coupling constant. It is given as 

Air 



in the leading order (LO), and it is 

■ NLO (Q 2 ) 



/3 ln(Q 2 /A 2 ) 



(A.l) 



Air 



a; 



/3 m(Q 2 /A 2 ) 



fa lnln(g 2 /A 2 ; 
/3 2 ln(Q 2 /A 2 ) 



(A.2) 



in the next-to-leading-order (NLO) 0. The renormalization scheme is MS. A is the 
QCD scale parameter which depends on the renormalization scheme . (3q and f3\ are 
the expansion coefficients of the /^-function and they are given by 



—C G - —TftNf , 
3 3 H f ' 



Pi 



—C 2 G - —C G N f 
3 3 1 



2C F N f , 

where C G , Cp, and T R are constants associated with the color SU(3) group: 

A^ 2 - 1 „ 1 



(A3) 



a 



G 



N r . 



2N C " 2 

A^c is the number of color (A r c =3) and Nf is the number of flavor. 



(A4) 



Appendix B. Splitting function 

and coefficient function 



Spin-independent flavor-nonsinglet splitting functions and the coefficient functions 
for calculating Fi, F2, and F3 are given in the following. The LO nonsinglet splitting 
function is given by 



P$l(x) = C F 



(B.l) 



where the integral including the function 1/(1 — x) + is defined by 



dx 



(l-x). 



dx 



1 — X 



(B.2) 
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The NLO nonsinglet splitting function for a ll q+q type" distribution [q NS = a.i(qi + qi)} 

i 

is given by Herrod and Wada M 



CU P F (x) + P A (x) + 5(1 - x) 



+ ~C F C A {P G (x)-P A (x) + 5(l-x) 
+ C F T R N f {p NF (x)-6(l-x)^ + ^} 



^-^ 2 + C(3)-85(oo) 

Y 2 + yvr 2 - C(3) + 8S(oc) 



(B-3) 



where P F (x), Pg(x), P Nf (x), and Pa(x) are given by Curci, Furmanski, and Petronzio 
ias 



Pf(x) 
Pg{x) 
P Nf {x) 
P A (x) 



1 + x 2 , , . . 

-2- lnxmfl — x) 

1 — x 

1+x 2 r 



1 — X 



+ 2x 



(1-x). 



1 + 



(1-x). 



,2 11 , 67 1 2 

In xH mxH ir 

3 9 3 



Y-lnx- -2(1 -x) 



lnx — -(1 + x) In 2 x —5(1 — x) 
40 

+ 2(1 + x)lnx + — (1 -x) , 



1 + x 2 dz , 1 - z 



In 



OO 



+ 2(1 +x) lnx + 4(1 -x) 



(B-4) 



1+X Jx/(l+x) z 

The C function is defined by ((k) = — and the numerical value (£(3)=1. 2020569. 



n=l 



rr 



is taken from Ref. 5(00) is given by the ( function as S(oo) = — |C(3) in the 

Gonzalez- Arroyo paper ||. In the case of a u q — q type", we simply change the sign of 
P A (P A ^-P A )[§. 



Pns-( x ) 



C 2 F { P F (x) - P A (x) + 6(1 - x) 



~ " + C(3) - 85(oo) 



+ -C F C A |P G (x)+P A (x) + 5(l-x) 
+ C7 J ,T fl JV / |p JVF (ar)-<5(l-a;)(i + ^7r 



^ + yvr 2 - C(3) + 8S(oc) 



(B-5) 



The coefficient functions are needed in calculating structure functions with the NLO 
corrections. The coefficient functions B N (x) for the structure functions, Pi, P2/X, and 
P 3 , are 



Pf s (x) 
P 2 ^(x) 



^[^)-4x] , 

C F F q (x) , 
C F [P 9 (x) - 2 -2s] 
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(B.6) 



where the function F q (x) is given by 



F q (x) 



3 1 + x2 l/o r X n l + I 2 , 2 , Tln(l-X) 

■-- + (9+5x)-2- lnx+2 1+x 2 -J ^ 

2(1 — x U 2 1 — x 1 — x 



-5(l-x)(9+-vr 2 ) 
(5.7) 



Appendix C. Laguerre polynomials 

and the LO evolution operator 



We use the Laguerre polynomials defined by [16 



n\x' m 



\ m ml 



m=0 



It should be noted that a slightly different definition, the above equation multiplied by 
n\, is sometimes used. A function of the variable x is first expressed as the function of 
x' = — In x, then it is expanded by the Laguerre polynomials. 

/(i = e" V ) = E/nA»M , (C.2) 

n 

where f n is the expansion coefficient given by the integral 

f n = I™ dx'e- x 'L n (x')f(e- x ') . (C3) 
Jo 

This integral is calculated numerically by a Gauss quadrature. If the initial distribution 
can be expressed in a simple analytical form f(x) = ax b (l — x) c , the coefficients are 
obtained by || 

with the recurrence relations go(c) = 1, g>jt+i( c ) — ~9k( c ){ c — k)/(k+ 1). 

The LO evolution operator is expanded by the Laguerre polynomials, and the 
method of calculating the expansion coefficients are discussed in Refs. j|, |], 0- We 
only show the results in the following. The coefficient is given by 

N j.k 
fe=0 K! 
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where B^ is calculated by the recurrence relation 

m— 1 

Bt +1 = Y.( P ™^-Prn-i-l)B* , (C.6) 

with 

5? = 1, B}=j2(Pj-Pj-i), B* = B* = ... = B k k _ 1 = . (C.7) 

3=1 

Namely, we first calculate the expansion coefficients of the splitting function (P n ), then 
B^ is calculated by the recurrence relation in Eqs. (C.6) and (C.7). 
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Figure Captions 



Fig. 1 Valence-quark distributions calculated by our program LAG2NS. The valence- 
quark distribution at Q 2 =4 GeV 2 is the HMRS-B distribution, which is evolved 
to the ones at Q 2 =200 GeV 2 . Three results are shown; the dotted, dashed, and 
solid curves are obtained by taking ten, twenty, and thirty Laguerre polynomials 
respectively. See text for the details. 

Fig. 2 Our evolution results are compared with CDHSW experimental data. We 
choose the CCFR distribution at Q 2 =3 GeV 2 in both LO and NLO cases, then we 
calculated Q 2 variations. F 3 structure functions at a;=0.045, 0.225, and 0.55 are 
shown. The solid (dashed) curves are the results of the NLO (LO) Q 2 evolution. 
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TEST RUN OUTPUT 



****** ENTER ****** F0,IORDER,ITYPE,IMORP,Q2,Q02,DLAM, 

NMAX,IQN,ILOG,XMIN,XMAX,NXSTEP 
F0= 4.0 IORDER= 2 ITYPE= 4 IMORP= 1 Q02= 4.0 Q2=200.000 DLAM= 0.190 
NMAX=30 IQN=1 ILOG=2 XMIN= 0.0100 XMAX= 1.00 NXSTEP= 30 
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This figure "figl-l.png" is available in "png" format from: 



http://arXiv.org/ps/hep-ph/9409289vl 



This figure "figl-2.png" is available in "png" format from: 



http://arXiv.org/ps/hep-ph/9409289vl 



